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Abstract.  Probabilistic  aspects  of  regional  ocean  model  pre¬ 
dictability  is  analyzed  using  the  probability  density  function 
(PDF)  of  the  irreversible  predictability  time  (IPT)  (called  r- 
PDF)  computed  from  an  unconstrained  ensemble  of  stochas¬ 
tic  perturbations  in  initial  conditions,  winds,  and  open 
boundary  conditions.  Two- attractors  (a  chaotic  attractor  and 
a  small-amplitude  stable  limit  cycle)  are  found  in  the  wind- 
driven  circulation.  Relationship  between  attractor’s  resi¬ 
dence  time  and  IPT  determines  the  t-PDF  for  the  short  (up 
to  several  weeks)  and  intermediate  (up  to  two  months)  pre¬ 
dictions.  The  r-PDF  is  usually  non-Gaussian  but  not  multi¬ 
modal  for  red-noise  perturbations  in  initial  conditions  and 
perturbations  in  the  wind  and  open  boundary  conditions.  Bi¬ 
furcation  of  r-PDF  occurs  as  the  tolerance  level  varies.  Gen¬ 
erally,  extremely  successful  predictions  (corresponding  to 
the  r -PDF’s  tail  toward  large  IPT  domain)  are  not  outliers 
and  share  the  same  statistics  as  a  whole  ensemble  of  predic¬ 
tions. 


1  Introduction 

Estimate  of  predictability  skill  of  regional  ocean  models  is 
important  but  difficult.  The  regional  ocean  dynamics  is  sen¬ 
sitive  to  variability  of  external  forcing  such  as  winds,  fresh 
water  discharge,  buoyancy  flux,  and  to  multi-scale  interac¬ 
tions  between  shelf  and  abyssal  currents  (Robinson  et  al., 
1996;  Lozano  et  al.,  1996;  Jiang  and  Malanotte-Rizzoli, 
1999;  Chu  et  al.,  1999a,  b,  c,  d,  2001;  Auclair  et  al.,  2003 
among  others).  Bathymetry,  coastlines,  and  physical  pro¬ 
cesses  such  as  wind  bursts,  fresh  water  discharge,  tides, 
storm  surges,  shelf  waves,  and  nonlinear  processes  affect 
coastal  and  abyssal  currents  and  may  cause  multi-attractor 
circulation  in  marginal  seas  where  robust  circulation  regimes 
(attractors)  are  more  complicated  than  simply  periodic  or 
even  quasi-periodic  ones. 
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Multiple  attractors  in  atmospheric  models  are  used  by  sev¬ 
eral  investigators  to  understand  the  dynamics  of  planetary 
flow  regimes  (Charney  and  De  Vore,  1979;  Molteni,  1996). 
Multi-modality  is  found  in  wind-driven  ocean  circulation 
models  (Schmeits  and  Dijkstra,  2000)  such  as  the  Kuroshio 
path  variation  south  of  Japan  (Masuda  et  al.,  1999).  The 
multi- modality  in  the  probability  density  function  (PDF)  of 
prediction  error  is  caused  by  the  spatial  phase  organization 
of  the  local  error  growth  rate  (Miller  and  Ehret,  2002).  Usu¬ 
ally,  the  multi-modality  represents  the  intermittency  of  the 
predictability  of  atmospheric  and  ocean  models  and  the  for¬ 
mation  of  non-Gaussian  PDF  of  the  prediction  error. 

Even  for  prediction  error  with  initial  Gaussian  PDF,  the 
PDF  of  predictability  time  scale  is  often  non-Gaussian  with  a 
tail  stretching  into  large  predictability  time  scale  values  (Chu 
and  Ivanov,  2002,  2004a,  b;  Chu  et  al.,  2002a,  b;  Denholm- 
Price,  2003).  Rare  forecasts  of  enhanced  predictability  dom¬ 
inate  these  tails.  Obviously,  if  the  occurrences  of  these  pre¬ 
dictions  increase,  the  ensemble  prediction  will  be  improved. 

In  this  study,  a  wind-driven  circulation  model  based  on 
the  barotropic  version  of  Princeton  Ocean  Model  (POM) 
(Blumberg  and  Mellor,  1987)  is  implemented  for  a  semi- 
enclosed  rectangular  basin  with  flat  bottom  to  reproduce  the 
multi-attractor  regime  and  to  study  the  variability  of  model 
predictability- skill  with  uncertain  initial  conditions,  external 
forcing  (wind)  and  open  boundary  conditions.  A  recently 
proposed  irreversible  predictability  time  (IPT)  r  (Ivanov  et 
al.,  1994;  Chu  et  al.,  2002a,  b)  is  taken  as  a  quantitative  mea¬ 
sure  of  model  predictability  skill. 

There  are  three  major  tasks  in  this  study:  (1)  to  find  mech¬ 
anisms  for  the  formation  of  non-Gaussian  PDF  of  IPT  (called 
r-PDF)  with  a  tail  stretching  into  large  IPT  domain,  (2)  to  de¬ 
tect  the  effect  of  stochastic  perturbations  in  external  forcing 
on  the  model  predictability,  and  (3)  to  identify  the  control 
parameter  (or  control  parameters)  for  the  variability  or  bifur¬ 
cation  of  the  model  predictability  skill  such  as  the  change 
of  symmetric  r-PDF  into  asymmetric  r-PDF  with  a  long 
tail  toward  large  IPT  domain.  With  fulfilling  these  tasks, 
the  following  fundamental  questions  will  be  answered:  Do 
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Fig.  1.  Two  areas  for  the  POM  integration:  (a)  domain- A  with  three 
rigid  boundaries  and  one  open  boundary  T' ,  and  (b)  domain-B  with 
four  rigid  boundaries. 


extremely  successful  predictions  (long  IPT)  share  the  same 
statistical  properties  as  predictions  of  short  and  intermediate 
durations?  How  do  the  stochastic  perturbations  in  external 
forcing  switch  the  attractors  (noise-induced  escapes)?  How 
do  they  affect  a  multi-attractor  system?  What  is  the  most 
important  parameter  to  determine  the  r-PDF  characteristics? 

The  paper  is  organized  as  follows.  Section  2  describes 
the  reference  solution  presenting  two  circulation  attractors: 
a  small-amplitude  stable  limit  cycle  and  a  chaotic  attractor. 
Section  3  discusses  the  method  for  generating  unconstrained 
ensembles  of  perturbations  added  to  the  initial  conditions, 
wind  forcing,  and  open  boundary  conditions.  Sections  4  and 
5  discuss  the  first  and  second  kinds  of  predictability.  Sec¬ 
tion  6  presents  the  conclusions. 


2  The  reference  solution 

2. 1  Numerical  modeling 

2.1.1  Two  domains 

Consider  a  flat  bay  centered  at  35°  N  and  bounded  by  three 
rigid  boundaries  (Domain- A  in  Fig.  la).  This  bay  expands 
1000  km  (1050  km)  in  the  north-south  (east-west)  direction. 
The  northern,  southern,  and  western  boundaries  are  rigid, 
and  the  eastern  boundary  is  open.  The  Cartesian  coordinate 
system  is  chosen  with  the  origin  at  the  southwest  comer.  The 
x-  and  y-axes  point  towards  the  east  and  north,  respectively 
(Fig.  1).  The  eastern  boundary  of  Domain-B  is  connected  to  a 
domain  with  three  rigid  boundaries  forming  a  closed  rectan¬ 
gular  domain  (Fig.  lb),  called  Domain-B  (Chu  et  al.,  1997). 

The  circulation  in  the  bay  is  modeled  with  the  barotropic 
version  of  POM.  The  vertical  mixing  and  vertical  viscosity 
are  not  included.  The  horizontal  mixing  is  parameterized  by 


the  bi-harmonic  operator.  The  bottom  stress  is  parameterized 
using  the  quadratic  drag  relation.  The  POM  was  evaluated 
thoroughly  using  the  observational  data  (Chu  et  al.,  2001). 

2.1.2  Integration  over  the  closed  domain  (Domain-B) 


The  POM  model  was  integrated  for  Domain-B  with  four  rigid 
boundaries  from  rest  and  zero  surface  elevation  (f ), 

(m,  fi,f)  =  0,  (1) 


and  forced  by  the  zonal  surface  pseudo  wind  stress  varying 
with  latitude 


Pa 


(2) 


where  wo=!0~3  m2  s-2,  is  the  mean  pseudo  wind  stress,  and 
pa  =  1.29  kg  m-3,  is  the  air  density  at  the  ocean  surface.  The 
time  step  is  chosen  as  2  min.  The  horizontal  resolution  is 
50  km. 


2.1.3  Integration  over  the  open  domain  (Domain- A) 

We  integrate  the  POM  over  Domain-B  with  four  rigid  bound¬ 
aries  (known  boundary  conditions)  from  the  initial  condi¬ 
tions  (1)  and  surface  forcing  (2)  and  take  the  solution  along 
the  middle  of  Domain-B  (v=1050  km) 

® 

as  the  reference  open  boundary  condition  for  the  Domain-A 
integration.  The  velocity  at  v= 1050  km  for  the  Domain-B  run 
is  nearly  zonal.  The  solution  for  the  Domain-B  integration  at 
day- 10  are  taken  as  the  initial  and  open  boundary  conditions 
for  the  Domain-A  integration, 

(u,  V,  f)  =  (u,  v,f)\t=\0day  (4) 

The  model  for  Domain-A  is  integrated  with  the  wind  forcing 
(2),  open  boundary  condition  (3),  and  from  the  initial  condi¬ 
tion  (4)  for  another  750  days.  The  initial  condition  is  an  un¬ 
closed  single  gyre  with  velocities  are  under  0.35 ms-1  and 
sea  surface  elevation  around  0.05  to  0.1  m.  Figure  2b  shows 
the  normal  velocity  along  the  open  boundary. 

2.2  Two-attractor  circulation 


Two  different  attractors  in  the  numerical  solutions  with  the 
same  external  forcing  and  dissipation  are  detected  by  the 
temporal  variability  of  total  kinetic  energy,  Lyapunov  dimen¬ 
sionality,  and  local  Lyapunov  exponents  (Anishenko,  1997). 

An  attractor  can  be  detected  either  in  the  state  space 
(Berloff  and  McWilliams,  1999)  or  in  the  phase  space.  As 
indicated  in  Appendix- A,  the  compressible  flow  ( u ,  v)  with 
small  sea  surface  elevation  can  be  represented  approximately 
by  the  geostrophic  stream  function  (see  Eq.  A2), 

K 

y{x,  y,  t)  =  4 >h(x,  y )  +  L  Ak(t)'i)k(x,  y ),  (5) 

k= 1 
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(a)  (b) 


Fig.  3.  Typical  IPT  landscape  projected  onto  error  phase  space 
|  Zp  | :  (a)  chaotic  attractor  and  (b)  limit  cycle. 


Fig.  2.  Initial  condition  for  modeling  (a)  and  open  boundary  condi¬ 
tion  U\j. 


where  is  the  harmonic  function  explained  in  Appendix 
A. 

The  orthonormal  functions  {4^}  (modes)  generate  a  K- 
dimension  phase  space.  The  series  (5)  is  truncated  at  ^=100. 
An  attractor  is  characterized  by  its  basin  volume  and  stabil¬ 
ity. 

2.2.1  Basin  volume 


If  the  model  dynamics  is  high-dimensional  and  dissipative,  it 
is  expected  that  generally  the  attractor  basin  is  not  a  smoothly 
connected  object,  but  a  fractal  or  a  riddled  basin  (Guchen- 
heimer  and  Holmes,  1983).  Accordingly  to  Kaneko  (1998), 
the  basin  volume  of  an  attractor  is  estimated  as  the  portion 
of  the  randomly  selected  initial  points  that  are  attracted  to 
it.  After  model  integration  over  100  days  with  10000  ini¬ 
tial  conditions,  the  attractor  is  identified  for  the  model  tra¬ 
jectory  falls.  Perturbed  model  trajectories  are  approximately 
distributed  between  two  attractors  in  the  proportion  as  1/9. 

A  manifold  of  initial  conditions  A^ito)  corresponding  to 
one  of  the  attractors  with  the  basin  volume  occupied  approx¬ 
imately  10%  of  phase  space  (the  attractor  is  identified  below 
as  a  small  amplitude  limit  cycle)  can  be  specified  through 
smoothing  (filtrating)  initial  conditions  accordingly  to 


A2 

Ak(to)  =  — r  2  Afc(fo), 

where  the  function  Xk  is  defined  by 


(6) 


Xk  =  Xomax  |A*(f0)| /Ak(t0), 

and  xo  is  the  filter  parameter;  is  the  eigenvalue  corre¬ 
sponding  to  the  basis  function  4^  (x,  y).  Chu  et  al.  (2003) 
demonstrate  that  the  filter  (6)  removes  high-order  modes  and 
low-order  small- amplitude  modes. 


2.2.2  Measure  of  attractor  stability 


To  estimate  attractor  stability  in  forecast  metrics 


J(f)  = 


ZTZ 


2 


(7) 


Z(x,  t)  =  4>(x,  t)  —  4>(x,  t),  Z(x,  to)  =  Zq,  (8) 


where  (x,  t)  denotes  the  spatial  and  temporal  coordinates; 
*(x,  t)  and  4>(x,  t)  are  the  reference  solution  and  an  indi¬ 
vidual  prediction,  respectively;  the  superscript  ‘T’  denotes 
the  transpose  operator,  ||  ||  is  the  Euclidian  norm,  we  use  the 
irreversible  predictability  time  determined  as 


T 


=  min 


J(t )  >  £2^  , 


(9) 


where  £=£o/  H'F  ||2,eo  is  tolerance  level  (accepted  prediction 
accuracy). 

IPT  differs  from  the  e-folding  time  (je)  or  doubling  time  if 
prediction  error  is  oscillatory  or  stochastic.  For  example,  the 
mean  IPT  over  an  ensemble  of  initial  conditions  is  defined 
by 


(r) 


Jit)  >£2)). 


(10) 


However,  the  e-folding  time  is  determined  by 


(re>  =  min  (f  |(  J(t))  >  e2)  . 


(ID 


with  the  corresponding  £2. 

Clearly,  unique  correlations  exist  among  the  landscape  of 
IPT  (dependence  of  IPT  on  magnitude  and  orientation  of  the 
initial  error),  the  local  Lyapunov  exponents  and  the  attractor 
types.  In  particular,  Chu  and  Ivanov  (2004a,  b)  pointed  out 
that  for  a  given  IPT  landscape,  the  signs  of  the  local  Lya¬ 
punov  exponents  can  be  determined. 

For  a  chaotic  attractor  with  at  least  one  positive  exponent 
and  one  negative  exponent,  the  IPT  landscape  has  a  crest-like 
structure  (Fig.  3a).  If  the  initial  error  is  oriented  perpendic¬ 
ular  to  the  crest,  the  prediction  error  grows  exponentially.  If 
the  initial  error  is  oriented  along  the  crest,  the  most  unstable 
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Fig.  4.  The  chaotic  attractor:  (a)  mean  kinetic  energy  normalized  by 
the  value  of  maximum  mean  kinetic  energy  of  limit  cycle  Eq  ;  (b) 
snapshot  of  stream  function  corresponding  to  a  maximum  energy 
state  of  attractor. 


Fig.  5.  The  small- amplitude  limit  cycle:  (a)  mean  kinetic  energy 
normalized  by  its  maximum  value  Eq ,  (b)  snapshot  of  stream  func¬ 
tion  averaged  over  time. 


Lyapunov  vector  decays  (stagnates).  A  hill-like  landscape 
corresponds  to  a  stable  limit  cycle  (Fig.  3b).  Here,  there  is 
no  specified  orientation  in  the  phase  space. 

2 . 3  Attractor  patterns 

2.3.1  Chaotic  attractor 

After  integrating  the  model  from  the  original  (taken  as  ac¬ 
curate)  initial  condition  (Fig.  2b),  the  attractor  is  character¬ 
ized  by  a-periodic  behavior  of  kinetic  energy  (Fig.  4a).  The 
streamfunction  with  the  maximum  energy  (Fig.  4b)  shows  a 
strong  cyclonic  basin-scale  gyre  with  size  of  about  350  km 
(left  upper  corner)  and  speed  of  about  0.9-1 .0ms-1.  Several 
smaller  size  (secondary)  eddies  occur  near  the  rigid  boundary 
and  in  the  center  of  the  basin.  The  highest  surface  elevation 
reaches  1  m. 

The  nonlinear  regime  is  identified  in  the  basin  by  two  non- 
dimensional  Rossby  numbers  (Pedlosky,  1987), 

U  1 

Ro i  =  max( - )  ~  0.1,  R02  =  max( - )  ~  0.12  (12) 

fL  fT 

where  /~  10-4  s-1  is  the  Coriolis  parameter;  t/~  1  ms-1  is 
the  characteristic  velocity  in  the  basin;  L~100km  and  r~l 
day  are  the  characteristic  spatial  and  temporal  scales  of  the 
flow,  respectively. 

Temporal  variation  of  the  kinetic  energy  averaged  over 
the  whole  basin  (Fig.  4a)  shows  a-periodic  oscillations  and 
transitions  between  meta-stable  states.  These  oscillations  are 
caused  by  interactions  between  the  basin- scale  cyclonic  gyre 
and  the  secondary  eddies  which  are  generated  by  barotropic 
instability  of  the  gyre  when  the  circulation  velocity  exceeds 
0.7-0. 8  ms-1.  The  instability  produces  the  secondary  eddies 
disposing  along  the  rigid  boundary  and  in  the  center  of  the 
basin  and  periodically  dissipating  due  to  viscosity  damping. 
Note  that  for  zero  horizontal  viscosity  the  chaotic  attractor 


does  not  develop.  Instead  an  unstable  limit  cycle  occurs.  This 
case  will  be  discussed  in  Chu  and  Ivanov  (2004a,  b). 

The  analysis  of  IPT  landscape  (Chu  and  Ivanov,  2004a,  b) 
demonstrates  the  existence  of  at  least  one  positive  and  one 
negative  Lyapunov  exponent.  It  identifies  the  observed  dy¬ 
namical  regime  as  a  chaotic  attractor. 

2.3.2  Stable  limit  cycle 

If  the  nonlinear  filtration  procedure  (4)  is  applied  to  the  ini¬ 
tial  condition  (4),  an  oscillated  stable  cyclonic  gyre  quickly 
develops  for  the  semi-enclosed  basin.  Figure  5a  shows  the 
temporal  variation  of  the  kinetic  energy  averaged  over  the 
whole  basin.  The  oscillation  period  is  about  150  days.  The 
nondimensional  amplitude  of  the  oscillation  is  around  0.1. 
Figure  5b  shows  temporally  averaged  stream  function. 

Such  a  circulation  pattern  is  represented  as  a  stable  limit 
cycle  (small  amplitude)  in  the  phase  space.  This  limit  cy¬ 
cle  is  identified  using  the  method  proposed  by  Berloff  and 
Meacham  (1997). 

3  Unconstrained  ensemble  of  perturbations 

Since  the  barotropic  model  is  computationally  efficient,  the 
stochastic  stability  of  the  reference  solution  to  unconstrained 
perturbations  of  the  initial  condition,  open  boundary  con¬ 
dition,  and  wind  forcing  is  studied  using  the  Monte-Carlo 
method  with  1000  ensemble  samples. 

Atmospheric  ensemble  modeling  shows  that  for  small 
samples  the  initial  error  can  not  cover  the  phase  space  homo¬ 
geneously.  It  leads  to  an  unrealistically-weak  growth  of  pre¬ 
diction  error  variance.  Therefore,  for  short  period  prediction 
it  requires  the  best  possible  knowledge  of  the  initial  state  in 
the  sub-region  of  the  phase  space  where  the  prediction  errors 
are  most  likely  to  grow.  In  different  forecast  centers,  such  as 
the  European  Centre  for  Medium-Range  Weather  Forecasts 
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and  the  National  Centers  for  Environmental  Prediction,  the 
directions  corresponding  to  the  maximum  error  growth  in  the 
phase  space  are  identified  by  the  singular  vectors  (Buizza  et 
al.,  1999)  and  the  breed  vectors  (Toth  and  Kalnay,  1997). 

This  study  utilizes  a  simple  Monte-Carlo  technique  for  the 
three  reasons.  First,  little  difference  in  the  statistics  is  found 
among  ensembles  of  1000,  5000,  10000  and  20000  sam¬ 
ples.  The  numerical  experiments  with  1000  ensemble  sam¬ 
ple  are  reliable.  Second,  a  Latin  hypercube  sampling  pro¬ 
cedure  (LHST,  2001)  is  used  to  increase  the  degree  of  uni¬ 
form  coverage  of  the  phase  space  by  initial  errors.  Third,  a 
non-parametrical  technique  using  the  Epanechinikov  kernel 
is  used  to  reconstruct  PDF  from  histogram.  Details  of  such  a 
reconstruction  can  be  found  in  Good  (1996). 

3.1  Perturbations  to  initial  condition 


Sensitivity  of  the  reference  solution  to  Gaussian-type  pertur¬ 
bations  is  commonly  studied  in  the  atmospheric  and  oceanic 
applications  (Robinson  et  al.,  1996).  In  this  study,  2D 
stochastic  perturbations  with  the  correlation  function 


G(r  —  r')  =  /2exp 


(r-rQ2 

R 2 


(13) 


are  used.  Here,  r=vi+yj,  is  the  displacement  vector;  (i,  j) 
are  unit  vectors  along  x  and  y  axes;  (R,  7)  are  the  correlation 
radius  and  intensity  of  the  noise.  The  technical  details  of 
generating  these  noises  can  be  found  in  Sabel’feld  (1991). 


3.2  Perturbations  to  wind  forcing 


Surface  winds  are  perturbed  by 

u'  =  [ u'(x ,  y,  t),  v'(x,  y,  t)]  =  [mi  (t),  M2 (t)]g(x,  y ),  (14) 


where  /jL\  and  fi2  are  independent  white  noise  processes  with 
zero  mean  and  variance  of  a2.  In  numerical  integration,  the 
noises  (/xi,/X2)  are  introduced  every  hour.  The  spatial  struc¬ 
ture  function  of  errors,  g,  is  parameterized  by, 


g(x,y)  =a 


7tXxXyerf 


exp 


/(x-x0)2  (y-yo)2) 

\  %  q  ) 


erf 


(15) 


which  describes  the  impact  of  the  localized  atmospheric  eddy 
activity  near  (xq,  yo)  on  the  surface  wind  perturbations  (Sura 
et  al.,  2001).  Here,  erf  is  the  error  function;  a  is  a  scaling 
parameter;  Xy)  are  the  de-correlation  scales.  Noise  in 
surface  winds  with  a 2 =28  m2  s-2  corresponds  to  observed 
typical  atmosphere  conditions  in  the  North  Atlantic  region 
(Wright,  1988). 


3.3  Perturbations  to  open  boundary  condition 


The  normal  velocity  of  the  reference  solution  along  the  open 
boundary  (ub)  is  perturbed  by 

u'b  =  B(t)  sm(jry/Ly),  (16) 


Fig.  6.  r-PDF  for  the  predictability  of  the  first  kind  with  white 
initial  noise  (curve-2),  and  red  noise  with  the  correlation  radius 
of  14  km  (curve-1)  and  70  km  (curve-3).  Here,  I2- 0.01  m2/s2  and 


where  B{t)  is  white  noise  with  zero  mean  and  variance  of 
82.  The  velocity  perturbation  ( u'b )  affects  inflow  and  outflow 
across  the  open  boundary. 


4  Predictability  of  the  first  kind 


Keeping  wind  forcing  (2)  and  open  boundary  condition  (4) 
accurate,  and  perturbing  initial  condition  (3)  with  white  or 
red  noise  represented  by  two  parameters  (variance  of  initial 
error  72  and  correlation  radius  R ),  the  predictability  of  the 
first  kind  (Lorenz,  1982;  Chu,  1999)  can  be  investigated  with 
various  combinations  of  (72,  R ,  £(2).  In  the  phase  space,  per¬ 
turbed  model  trajectories  start  from  the  basin  of  chaotic  at¬ 
tractor,  i.e.  from  the  original  initial  condition  (3).  The  r-PDF 
is  sensitive  to  the  variation  of  £q. 

For  a  given  tolerance  level  £q  (0.5),  the  r-PDF  is  close  to 
Gaussian  shape  (skewness  of  0.08)  for  the  white  noise  (R= 0, 
curve- 1  in  Fig.  6).  Clearly,  any  model  trajectory  starting  from 
the  attractor  basin  does  not  leave  the  attractor  and  the  IPT 
statistics  is  only  determined  by  the  stability  of  this  chaotic 
attractor. 

When  the  correlation  radius  R  increases,  the  volume  in 
the  phase  space  covered  by  the  initial  error  increases  (Chu 
and  Ivanov,  2004a),  and  the  two-attractor  structure  of  the 
reference  circulation  occurs.  For  a  given  error  intensity, 
72=0.01  m2  s-2,  the  r-PDF  is  near  Gaussian  (curve-1  in 
Fig.  6,  skewness  near  0)  for  R= 0  km,  non-Gaussian  (curves-2 
in  Fig.  6,  skewness  of  0.58)  for  R=14km,  and  non-Gaussian 
with  a  tail  stretching  to  large  IPT  domain  (curves-3  in  Fig.  6, 
skewness  of  0.80)  for  7?=70km. 
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Fig.  7.  r-PDF  for  the  predictability  of  the  second  kind  (wind  forc¬ 
ing)  with  a  given  wind  error  intensity  (cr2=  42  m2  s-2)  and  various 
tolerance  levels  (£q):  (a)  small  tolerance  levels:  0.005  (curve- 1)  , 
0.01  (curve-2),  and  (b)  intermediate  and  large  tolerance  levels:  0.05 
(curve-1),  0.1  (curve-2)  and  0.5  (curve-3). 


5  Predictability  of  the  second  kind 

Keep  the  initial  condition  (3)  accurate,  and  perturb  the  wind 
forcing  (2)  and  open  boundary  condition  (4)  with  white 
noises.  The  noise  the  intensity  is  a2  for  the  wind  forcing  and 
8 2  for  the  open  boundary  condition.  The  second  kind  of  pre¬ 
dictability  (Lorenz,  1982)  can  be  investigated  with  various 
combinations  of  (a2,  sfy  and  (82,  £q).  Since  the  predictabil¬ 
ity  of  high-resolution  regional  models  depends  crucially  on 
the  wind  forcing  (Robinson  et  al.,  1996;  Chu  et  al.,  1999a, 
b,  c;  Burillo  et  al.,  2002  and  others)  and  the  specified  nor¬ 
mal  velocity  at  the  open  boundaries,  and  since  errors  in  the 
open  boundary  condition  cause  a  rapid  decrease  of  model 
predictability  skill  (Chu,  1999;  Jiang  and  Malanotte-Rizzoli, 
1999),  investigation  of  the  second  kind  of  predictability  is  of 
great  importance  for  regional  ocean  modeling. 

5.1  Residence  time  and  r-PDF 

In  contrast  to  the  predictability  of  the  first  kind,  analysis  of 
the  second  kind  of  predictability  is  more  difficult  because  of 


the  multi-attractor  system.  “Noise  induced  escape”  occurs 
when  a  model  trajectory  switches  over  attractors  (Anishenko, 
1997).  Therefore,  in  addition  to  IPT,  the  residence  time  ( xres ) 
should  be  used  with  the  definition  of  a  maximum  time  period 
for  a  perturbed  trajectory  keeping  in  the  same  attractor. 

In  general,  the  residence  time  depends  on  the  basin  volume 
of  the  attractor,  initial  position  of  model  trajectory  at  the  at¬ 
tractor,  dynamical  attractor  stability  and  noise  intensity  (a2 
or  82)  (Kaneko,  1997).  Since  there  is  no  common  method  to 
estimate  xres  on  the  attractor  for  a  high-dimensional  dissipa¬ 
tive  dynamical  system  (Kaneko,  1997),  a  practical  method  is 
proposed. 

For  a  small  £q,  it  is  expected  that 

T  <  Tres , 

which  shows  that  the  model  losses  the  predictability  skill  be¬ 
fore  the  model  trajectory  leaves  the  occupied  attractor.  The 
r-PDF  is  determined  by  the  stability  of  a  single  attractor  and 
should  be  near  Gaussian  type.  For  a  sufficiently  large  £q,  it 
is  expected  that 

r  >  zres , 

which  shows  that  the  model  keeps  the  predictability  skill  af¬ 
ter  the  model  trajectory  leaves  the  occupying  attractor.  This 
leads  to  the  change  of  the  r-PDF  structure.  Therefore,  switch 
among  attractors  can  be  determined  through  analyzing  the 
change  of  r-PDF  structure  with  increasing  £q. 

The  proposed  methodology,  of  course,  can  not  deter¬ 
mine  the  residence  time  xres  and  the  number  for  a  trajec¬ 
tory  switching  among  attractors.  However,  it  determines  the 
change  of  the  r-PDF  structure  that  links  to  the  duration  of 
the  perturbed  trajectory  along  different  attractors.  Besides,  it 
uniquely  identifies  the  trajectory  switching  from  the  chaotic 
attractor  to  the  limit  cycles  when  the  stochastic  perturbation 
is  added  to  external  forcing. 

5.2  Perturbations  to  wind  forcing 

The  two  parameters  o2  and  £q  determine  the  r-PDF  struc¬ 
ture.  For  small  tolerance  level  (£q<0.01)  and  strong  wind 
perturbations  (a2>40m2  s-2),  the  r-PDF  is  nearly  symmet¬ 
ric  with  large  variance  of  IPT.  The  kurtosis  and  skewness  are 
approximately  3.1  and  0.1  (Fig.  7a). 

With  increase  of  £q  (> 0.075),  the  r-PDF  becomes  asym¬ 
metric  with  smaller  IPT  variance  and  a  long  tail  stretching  to 
large  IPT  domain  (Fig.  7b).  The  skewness  is  0.60  for  curve- 1 
and  0.75  for  curve- 3.  This  shows  that  the  tolerance  level  £q  is 
a  controlling  parameter  for  model  predictability  skill.  The  r- 
PDF  is  symmetric  for  small  £q  and  asymmetric  for  large  £(2. 
With  a  relatively  large  tolerance  level  (£q=0.1),  the  asymme¬ 
try  of  the  r-PDF  structure  increases  (longer  tail  toward  large 
IPT  domain)  with  the  increase  of  noise  intensity  a2  (Fig.  8). 
The  extra- successful  predictions  (Chu  et  al.,  2002c)  corre¬ 
sponding  to  the  PDF  tail  (toward  the  large  IPT  domain)  are 
not  outliers  since  they  share  the  same  statistics  with  other 
predictions  with  short  and  intermediate  IPTs. 
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Fig.  8.  r-PDF  for  the  predictability  of  the  second  kind  (wind  forc¬ 
ing)  with  a  given  tolerance  level  £q=0.1  with  various  wind  error 
intensity:  o2- 14  m2  s-2  (curve-1),  28  m2  s-2  (curve-2),  42  m2  s-2 
(curve-3),  and  56  m2  s~2  (curve-4). 


Fig.  9.  r-PDF  for  the  second  kind  of  predictability  (open  bound¬ 
ary)  with  agiven  open  boundary  error  intensity  (<$2=0.2  m2  s-2)  and 
various  tolerance  levels  (Sq):  0.005  (curve-1),  0.01  (curve-2),  0.1 
(curve- 3)  and  0.5  (curve-4). 


Another  important  result  is  that  the  multi-attractor  struc¬ 
ture  of  the  modeled  circulation  does  not  necessary  lead  to  the 
multi-modal  r-PDF  since  the  perturbed  trajectory  can  switch 
among  different  attractors.  This  is  different  from  a  recent  re¬ 
sult  reported  by  Miller  and  Ehret  (2002)  that  if  the  initial  er¬ 
ror  co- variances  are  large  enough  such  that  a  given  perturbed 
trajectory  has  high  possibility  to  start  from  different  basin 
attractions.  The  initial  Gaussian  r-PDF  will  evolve  into  a 
bimodal  r-PDF. 

Existence  of  a  tail  stretching  into  large  IPT  domain  is  the 
most  common  feature  of  r-PDF  for  short  and  intermediate 
forecasts  (Fig.  8).  Increase  of  the  noise  intensity  a2  enhances 
the  asymmetry  of  the  r-PDF.  For  example,  the  skewness  is 
0.38  for  a2=14m2  s-2  (curve- 1)  and  0.96  for  a2=56m2  s-2 
(curve-4).  The  extremely  successful  predictions  (Chu  et  al., 
2002b)  corresponding  to  the  r-PDF  tails  are  not  outliers 
since  they  share  the  same  statistics  with  the  usual  predictions. 
This  will  be  pointed  out  in  Chu  et  al.  (2004b)  that  the  Weibull 
distribution  is  the  best  fit  for  the  r-PDFs  illustrated  in  Fig.  8. 

The  r-PDF  is  nearly  symmetric  for  small  £q  and  asym¬ 
metric  for  large  £q.  Loss  of  the  r-  PDF  symmetry  occurs  as 
£q  >0.075.  This  indicates  that  variation  of  the  tolerance  level 
£q  may  lead  to  a  drastic  change  of  statistics  for  measuring 
the  model  predictability  skill.  A  set  of  structural  change  of 
the  r-PDF  with  the  increase  of  £q  may  be  considered  as  the 
one-step  bifurcation  sequence: 

r-PDF  symmetry  — >  r-PDF  asymetry.  (17) 

5.3  Perturbations  to  open  boundary  condition 

For  small  tolerance  level  (£q<0.01),  the  r-PDF  is  already 
asymmetric  with  a  maximum  value  at  small  IPT  domain 
(about  5  days,  i.e.,  the  most  probable  IPT)  and  a  tail  stretch¬ 
ing  to  large  IPT  domain  up  to  40-50  days  (curves- 1  and 
curve-2  in  Fig.  9).  The  skewness  is  3.8  for  £2 =0.005  (curve- 
1)  and  3.9  for  £q  =0.01  (curve-2).  The  estimated  mean  res¬ 


idence  time  is  about  10  days,  the  perturbed  model  trajecto¬ 
ries  departure  quickly  from  the  chaotic  attractor  to  the  small 
amplitude  limit  cycle.  With  increase  of  £q  (>0.025),  the 
IPT  variance  and  the  degree  of  r-PDF  asymmetry  reduces 
(curves-3  and  curve-4  in  Fig.  9).  The  skewness  reduces  to 
0.7-0. 8.  This  is  caused  by  the  perturbed  trajectories  having 
time  to  return  to  the  chaotic  attractor  from  the  limit  cycle. 

For  some  values  of  (£q,  8 2),  there  exist  “outliers”  -  predic¬ 
tions  with  statistical  properties  differing  considerably  from 
those  of  the  prediction  ensemble.  For  example,  when  £q=0.  1, 
<52=0.3  m2  s-2,  the  r-PDF  shows  a  complex  structure  with 
two  distinct  features  for  r<40  day  and  r>40  day  (Fig.  10a). 
The  predictions  for  r<40  day  may  be  identified  as  outliers. 
Therefore  the  r-PDF  in  Fig.  10a  can  be  represented  as  the 
combination  of  two  r-PDFs:  (a)  uniform  distribution  (shown 
in  Fig.  10b)  for  12  h  <r<39  days  and  (b)  near  Gaussian  dis¬ 
tribution  with  the  mean  and  variance  equal  to  50  day  and  64 
day2,  respectively.  A  set  of  structural  change  of  the  r-PDF 
with  the  increase  of  £q  may  be  considered  as  the  two-step 
bifurcation  sequence: 

asymetric  r-PDF  with  large  variance  — > 

— >  symmetric  r-PDF  +  outliers  — > 

->  asymetric  r-PDF  with  small  variance.  (18) 

Existence  of  bifurcation  sequences  (17)  and  (18)  should  be 
checked  in  atmospheric  and  oceanic  models. 

6  Conclusions 

(1)  Several  predictability  regimes  in  two-attractor  circula¬ 
tions  (the  chaotic  attractor  and  small-amplitude  limit  cycle) 
reproduced  by  a  numerical  model  have  been  identified  in 
a  semi-closed  basin  with  flat  bottom.  The  r-PDF  is  non- 
Gaussian  for  intermediate  and  large  amplitude  prediction  er¬ 
rors,  and  near  Gaussian  for  small  amplitude  prediction  errors. 
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Fig.  10.  “Outliers”  -  forecasts  with  <$2=0.3m2  s-2,  £q=0.1:  (a)  r- 
histogram  computed  for  the  perturbed  open  boundary,  (b)  best  fit 
to  the  homogeneous  distribution  for  r<40  day,  and  (c)  best  fit  to 
Gaussian  distribution  for  t>40  day. 


(2)  For  the  predictability  of  the  first  kind,  the  error  evolu¬ 
tion  depends  strongly  on  the  correlation  radius  of  the  initial 
error  and  its  variance.  If  the  initial  co- variances  and  correla¬ 
tion  radius  are  large  enough  so  that  there  is  large  possibility 
that  a  given  perturbed  model  trajectory  does  not  start  from 
the  same  attractor  basin  as  the  mean  ensemble,  the  r-PDF  is 
asymmetric  (non-Gaussian)  with  a  tail  stretching  into  large 
IPT  domain. 

(3)  The  stochastic  perturbations  added  to  the  wind  forcing 
and  open  boundary  condition  lead  to  switch  of  the  perturbed 
model  trajectories  over  attractors  (i.e.,  the  noise-induced  es¬ 
cape).  It  is  found  that  the  perturbed  model  trajectories  switch 
from  the  chaotic  attractor  to  the  limit  cycles  and  switch 
back.  Relationship  between  the  residence  time  at  the  attrac¬ 
tor  and  IPT  determines  the  shape  of  r-PDF.  For  z<zres  per¬ 
turbed  model  trajectories  do  not  have  time  to  switch  from  the 
chaotic  attractor  to  the  limit  cycle.  This  causes  the  r-PDF 
near  Gaussian.  For  z~zres  the  perturbed  model  trajectories 
can  be  either  in  the  initial  chaotic  attractor  or  in  the  limit 
cycle.  This  causes  the  bimodal  (not  often)  or  asymmetric 
non-Gaussian  r-PDF.  For  z>zres ,  the  perturbed  trajectories 
have  sufficient  time  to  switch  to  the  limit  cycle,  and  cause 
asymmetric  r-PDF. 

(4)  The  perturbed  model  trajectories  attracted  to  the  limit 
cycle  share  the  same  statistical  properties  as  the  trajectories 
still  along  the  initial  chaotic  attractor.  In  general,  the  numer¬ 
ical  computations  do  not  indicate  the  outlier  feature  of  pre¬ 
dictions  except  for  some  forecasts  realized  for  specific  com¬ 
binations  of  the  tolerance  level  and  variance  of  perturbations 
added  to  the  open  boundary  condition. 

(5)  Two  different  bifurcation  sequences  in  changes  of  z- 
PDF  are  found  with  increase  of  £q.  For  uncertain  wind  forc¬ 
ing,  one-step  bifurcation  is  found  with  symmetric  r-PDF 
changing  to  asymmetric  r-PDF.  For  uncertain  open  bound¬ 
ary  conditions,  a  more  complicated  two-step  bifurcation  is 
found  in  the  model  predictability  skill. 


Appendix  A 


Spectral  decomposition 


A  two-dimensional  compressible  circulation  ( u ,  v)  can  be 
represented  by 


- T  — ,  v  =  +  — , 

dy  dx  dx  dy 


(Al) 


where  and  O  are  the  geostrophic  streamfunction  and  ve¬ 
locity  potential,  respectively.  Let  (T,  T')  be  the  rigid  and 
open  boundaries.  The  geostrophic  streamfunction  is  de¬ 
composed  into  (Eremeev  et  al.,  1992;  Chu  et  al.,  2003) 


00 

^(x,  y,  t )  =  VH(x,  y,t)  +  y2  A.k{t)^k{x,  y). 


k= 1 


(A2) 


where  {T^}  are  the  eigen-functions  of  the  spectral  problem  , 


A^k  (x,  y)  =  -Xk^k(x,  y),  %  |rur/  =  0,  (A3) 


with  Xk  the  eigen- values  of  the  horizontal  Laplacian  (A). 
The  component  is  the  harmonic  function  satisfying  the 
boundary  conditions 


lr  =  0  ,  ^ h  Ir  =  'P b , 


(A4) 


where  4^  is  computed  from  knowledge  of  the  normal  veloc¬ 
ity  at  open  boundary. 

The  velocity  potential  O  is  decomposed  into 

00 

<£>  =  y]Bm(r)Om(x,3;),  (A5) 

m= 1 

where  {Om}  are  the  eigen-functions  of  the  spectral  problem 
A(Fm  (jc,  y)  —  /xmOm(v,  y), 

d>m  lr'  =  0, 

n  •  VOm  |r  =0,  (A6) 

where  n  is  the  unit  vector  normal  to  Tr;  {/zm}  are  the  eigen¬ 
values. 

The  basis  functions  {4^,  Om}  generate  the  phase  space 
which  can  be  used  for  the  spectral  analysis  of  ocean  circu¬ 
lation  in  a  basin  with  real  coastlines.  Since  the  surface  eleva¬ 
tion  is  small  comparing  to  the  water  depth  H , 

S/H  «  1 

the  quasi-geostrophic  approximation  is  reasonable  for  the 
flow  reproduced  in  the  present  study  (Pedlosky,  1987).  Thus, 
in  the  main  text,  only  the  geostrophic  streamfunction  4*  is 
considered. 
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